Topology of Collisionless Relaxation 
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Using extensive molecular dynamics simulations we explore the fine-grained phase space structure 
of systems with long-range interactions. We find that if the initial phase space particle distribution 
has no holes, the final stationary distribution will also contain a compact simply connected region. 
The microscopic holes created by the filamentation of the initial distribution function are always 
restricted to the outer regions of the phase space. In general, for complex multi-level distributions 
it is very difficult to a priori predict the final stationary state without solving the full dynamical 
evolution. However, we show that for multi-level initial distributions satisfying a generalized virial 
condition, it is possible to predict the particle distribution in the final stationary state using Casimir 
invariants of the Vlasov dynamics. 
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Statistical mechanics of systems in which particles in- 
teract through long-range (LR) forces remains an out- 
standing challenge^!]. Unlike short-range systems, which 
relax to thermodynamic equilibrium through binary col- 
lisions, systems with LR interactions become trapped in 
the quasi-stationary states (qSS) the life time of which 
diverges with the number of particles (245|. The fun- 
damental difficulty in studying the qSS is that these 
states explicitly depend on the initial particle distribu- 
tion. Furthermore, the lack of ergodicity intrinsic to LR 
systems prevents application of the standard equilibrium 
statistical mechanics which has proven to be so powerful 
for describing the many-body systems with finite-range 
forces j6|. In spite of these difficulties, it has been recently 
observed that there is a significant degree of universality 
in the process of collisionless relaxation. Many different 
systems ranging from plasmas 7] to gravitational clus- 
ters [Sl, Isl have been found to relax to qSS with a char- 
acteristic core-halo structure. 

In thermodynamic limit, LR systems are collisionless 
- particles move under the action of the mean-field po- 
tential produced by all the other particles. In general, 
the mean-field potential has a complex dynamics, char- 
acterized by quasi-periodic oscillations. It is possible, 
therefore, for some particle to enter in resonance with the 
oscillations and to gain large amounts of energy at the 
expense of the collective motion. This process, known in 
plasma physics as the Landau damping [10|, diminishes 
the amplitude of oscillations and leads to the formation 
of a tenuous halo of highly energetic particles which sur- 
rounds the high density core. When all the oscillations 
die out a qSS will be born. The qSS is characterized by a 
broken ergodicity — since the only mechanism through 
which particles can gain or loose energy is the Landau 
damping, once the mean-field potential becomes station- 
ary, the dynamics of each particle becomes integrable (for 
spherically symmetric systems) [6|. When this happens 
there is no longer a mechanism through which highly 



energetic particles of the halo can equilibrate with the 
particles of the core. 

The relaxation to qSS is very similar to the process 
of evaporative cooling. As some particles enter in res- 
onance with the collective oscillations they gain energy, 
while cooling down the core region. The Hamiltonian 
dynamics of LR systems is governed by the collision- 
less Boltzmann (Vlasov) equation [llj . This equation re- 
quires that the one-particle distribution function evolves 
in the phase space as the density of an incompressible 
fluid. This means that the core region can not collapse 
to the minimum of the potential energy- since this would 
violate the incompressibility requirement imposed by the 
Vlasov flow - instead the maximum phase space den- 
sity can not exceed that of the initial distribution. For 
one-level initial distributions this observation allows us 
to predict the qSS without having to explicitly solve the 
Vlasov equation or perform the MD simulations. When 
all the oscillations have died out, the particles in the 
core should occupy all the low-energy states up to the 
maximum phase space density permitted by the original 
one-level waterbag distribution. The particle distribu- 
tion inside the core will then be the same as that of a 
fully degenerate Fermi gas [7|, |8|, [l^ ■ On the other hand 
the particles in the halo will be approximately uniformly 
distributed up to the maximum energy corresponding to 
the location of the parametric resonance. The Fermi en- 
ergy and the phase space density of the halo particles can 
then be obtained from the requirement of the conserva- 
tion of the norm and of the total energy of the system. 
The theory has been found to be extremely successful, 
allowing us to predict the distribution functions of con- 
fined plasmas [71 . Id and 2d gravitational systems J8|], 
the HMF model [l^I, etc. The question that we would 
like to address in this Letter is how to extend the theory 
described above to more complex initial particle distri- 
butions 13-151. 



To be specific we will study a class of distributions 



which are compact with a simply connected support (no- 
holes). The distributions have L different phase space 
density levels, see Fig. la. To demonstrate the theory we 
will use a paradigmatic model of a system with LR forces 
known as the Hamiltonian Mean-Field (HMF )model il,]. 
The HMF model describes N particles that arc con- 
strained to move on a circle of radius one. The dynamics 
is governed by the Hamiltonian 
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where the angle 9i is the position of i'th particle and 
Pi is its conjugate momentum [^, |l6|. The macroscopic 
behavior of the system is characterized by the magne- 
tization vector M = {Mx,My), where M^ = {cos 6), 
My = {sin 6), and (• • • ) stands for the average over all 
the particles. The Hamilton's equations of motion for 
each particle reduce to 



e, = -M^{t) sin 61, (t) + My{t) cos 6l.,(t). 



(2) 



Since the Hamiltonian does not have an explicit time 
dependence, the average energy per particle. 
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is conserved. For symmetric distributions {9 — s- —9) 
My = throughout the evolution, so that the macro- 
scopic dynamics is completely determined by M^(t), 
which for simplicity we will write as M(i) [12l.ll7[. 

Using the MD simulations, we first calculate the one- 
particle distribution function in the qSS, f{9,p). From 
Jean's theorem, in the stationary state f{9,p) — /(e), 
where e{9,p) = p^/2 -f 1 — Mcos9 is the one-particle 
energy and M = {cos 9) is the magnetization. The ini- 
tial phase space particle distribution is shown in Fig la. 
It consists of three different levels with the phase space 
densities - from inside to outside - iji , 772 , and 773 . To ob- 
tain the f{e) we run the simulation until the systems has 
reached a qSS. We then separate all the particles into 
bins of energy width de. To calculate the distribution 
function /(e) the fraction of the particles in each bin is 
divided by the density of states 
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S{s-p^/2-l + Mcos9)dpd9, (4) 



where S{x) is the Dirac delta function and the integral 
is performed over all the phase-space, —00 < p < -foo, 

-TT < 9 < IT. Using S[f{x)] = J:^H^ - X^)/\f'ix^)l 

where Xi are the roots of f{x), the integration in mo- 
mentum leads to 
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FIG. 1. (color online), (a) The initial 3 level phase space 
distribution. Colors are used to denote the difference phase 
space densities. There are a total of A'^ = 10^ particles dis- 
tributed over three three level with densities: rji — 0.39 (0 < 

bl < Pmax/S) in blue, 772 = 0.56 {pmax/3 < \p\ < 2pmax/3) 

in red, and r/a = 0.17 {2pmax/3 < \p\ < pmax) in green. 
The initial magnetization is Mo = 0.8 {9max ~ 1.131) and 
Pmax = 0.59. Panel (b) shows a snapshot of phase space at 
t = 200, demonstrating the mechanism of mixing through 
filamentation. Panel (c) shows a phase space snapshot at 
t — 20000, after the qSS has been estabUshed. Notice the 
characteristic core-halo structure of the particle distribution. 
In panel (d) we plot the total distribution function /(e) in the 
qSS (solid black curve) and the partial distribution functions 
/n(e) for each phase space level (dashed curves). 



where 9max — tt ii e > 1 + M, and 9„iax = cos~-'^[(l — 
e)/Af] if e < 1-1- M. Without loss of generality we will 
take A/ > 0. Performing the integral in Eq. ([5]) yields 
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if K < 1, 
if K > 1, 



(6) 



where k = (e — 1 -I- M)/2M, and K{k) is the complete 
elliptic integral of the first kind. Note that k < 1 (k > 1) 
corresponds to librating (rotating) orbits. The results of 
these calculations are shown in Fig. 1. The panel (a) of 
Fig. 1 shows the initial particle distribution. The stretch- 
ing and folding, characteristic of collisionless relaxation, 
appears as the filamentation of the original phase space, 
panel (b). Fig. Ic provides a snapshot of the phase space 
particle distribution in the qSS. Notice the appearance 
of the characteristic core-halo structure. In panel (d) 
we plot the total particle distribution function /(e) in 
the qSS, and the partial distribution functions /ri(e) for 
particles which at i = were inside the levels of different 
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FIG. 2. (color online). The volume fractions (j}n{£) of the 
phase space occupied by the 3 density levels in the final qSS: 
771 (blue long dashed curve), 772 (red dot-dashed curve), r/a 
(green short dashed curve). The black solid curve shows the 
volume fraction occupied by sum of all three density levels. 
Note that in the core region, which extends up to e ~ 0.4, all 
the phase-space is fully occupied. On the other hand the halo 
is 90% empty phase space. 



phase space densities ??„ . The distribution functions have 
a complicated structure which, unfortunately, shed very 
little light on the properties of the qSS. 

The Vlasov dynamics requires that the phase space 
area occupied by each density level of the initial distribu- 
tion function must be preserved throughout the dynami- 
cal evolution [18]. To explore this feature we next study 
the fraction of the phase space volume at energy e occu- 
pied by the level ??„, 0„(e) = /„(£)/??„ • In Fig. 2 we plot 
the 4>n{£) in the qSS for a 3- level distribution of Fig. 1. 
Once again we see a very complex mixing of the density 
levels over the phase space. Amazingly, however, when 
all the volume fractions are summed together, a very sim- 
ple core-halo structure emerges, the solid curve of Fig. 2. 
For energies up to e w 0.4, the phase space is completely 
occupied by the density levels of the original distribution 
function. There are no holes in the central " core" region. 
All the vacancies are confined to the outer "halo" region. 
Indeed, almost 90% of the halo is an empty phase space. 

The simplicity of the core-halo structure suggests that 
it might be possible to predict the final qSS a priori, 
without having to solve explicitly the full many-body dy- 
namics. This, indeed, was the case for the one-level wa- 
terbag distributions [1^ , for which the core was perfectly 
described by the fully degenerate Fermi-Dirac (FD) dis- 
tribution, while the halo energy could be calculated us- 
ing the theory of parametric resonances. Unfortunately 
the situation is much more complex for many-level sys- 
tems. To get a better feel for the dynamics leading to 
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FIG. 3. (color online). Schematic of phase space evolution 
for a virialized two-level distribution: panel (a) is the initial 
condition; panel (b) the fine-grained distribution function in 
the final qSS. Note that the central core region has no holes 
(white squares). For a virialized initial condition, the oscilla- 
tions are suppressed and the resonances are not excited. For 
such initial distributions the halo is populated only by the 
"outer" (green) level 



formation of qSS, in Fig. 3 we show a schematic evolu- 
tion of a two-level system from the initial to final state. 
The phase space is divided into macrocells of area dpd9. 
Each macrocell is then subdivided into microcells. The 
incompressibility of Vlasov dynamics requires that each 
microcell is occupied by at most one density level. In 
the qSS the central core region does not contain any mi- 
croscopic holes (white squares) which are all confined to 
the outside halo region. Although the evolution leads 
to a completely degenerate core (with no holes) it can 
not, in general, be described by the FD statistics. The 
problem is that the degenerate limit (T = 0) of the FD 
statistics requires that the low energy states must be oc- 
cupied by the levels with largest value of rj^. To be in the 
ground state a FD system must be stratified so that for 
e < £1 only levels rja should be present, for £1 < £ < £2 
there should be only levels ijh, for £2 < £ < £3 there 
should be only levels rjc etc., where rja > rji, > rjc > ... 
and £i,£2,£3... are the Fermi energies. This, however, is 
not the case, as can be clearly seen from Figs. 2 and 
3. There is a mixture of different levels inside the core 
region. Thus, the statistics of density levels can not, in 
general, be obtained a priori without explicitly solving 
the full many-body dynamics. 

The complex mixture of different density levels is a 
consequence of the parametric resonances produced by 
the particle-wave interactions. If the resonances can be 
suppressed, the structure of the qSS should be much sim- 
pler [6| . The resonances will not be excited if the initial 
particle distribution satisfies a generalized virial condi- 
tion (GVC). To derive the GVC we require that the mean 
square angle (0^) does not vary significantly with time, 
i.e. that there are no strong envelope oscillations. This 
will be the case if the two temporal derivatives of {d{t)'^) 
vanish, i.e. if (9p} = and (p^) — (cos 0) (6' sin 6') = 0. 



The GVC can be satisfied by adjusting the maximum 
momentum of the particles of each density level in the 
initial distribution function. 

In the absence of resonances there is no mechanism for 
the individual particles to gain energy. Therefore, the 
maximum one-particle energy in the qSS will be the same 
as the maximum one-particle energy in the initial distri- 
bution. For a L-level distribution function satisfying the 
GVC, the mixing will be restricted to the consecutive en- 
ergy levels [£i,£i+i], allowing us to write a simple ansatz 
solution for the distribution function: 

L 

4 = 1 

where {si} are the L + 1 threshold energies that sep- 
arate regions of different phase-space density. Xi is 
the the fraction of the phase space volume occupied by 
the level rji in the the phase-space region with energy 
[ei,£i+i]. We define 77^+1 ee and El+i = e^ax = 
p'^^^/2 + l — Mo cosOmax, which is the energy of the most 
energetic particle from the initial condition. Lack of res- 
onances permits the density level transfers only between 
the consecutive energy levels [£i,£i+i] and [ei+i,ei+2]- 
The conservation of the phase space volume occupied by 
each density level provides us with L coupled equations 

Vi + xiV2 = iyi/m, (8) 

(1 - Xt-i)V, + x^V,+i = vJt],, ^ = 2, L, (9) 

where Vi is the fraction of particles in the density level 
r]i, and Vt = J^' g{e) de is the total phase-space volume 
between the energies e^-i and e^. The minimum energy 
is £0 = 1 — Ms, where Ms is the magnetization of the 
qSS. 

The phase space volumes transfered from one energy 
level to the subsequent energy level must be conserved. 
This means that all the Xi's are related by Xt-i^i ~ 
XiVi+i. Finally, the conservation of the total energy per 
particle and the self-consistency requirement for magne- 
tization, 




FIG. 4. The position and the velocity distributions in the 
qSS state. The symbols are the results of MD simulations. 
The dashed red curves are the predictions of the LB theory, 
while the solid black curves are the prediction of the present 
theory. There are no adjustable parameters. The initial two- 
level distribution function satisfies the GVC. We consider two 
examples: panels (a) and (b) are the results for the initial 
distribution with ryi = 0.35 (0 < \p\ < Pmax/S), 7)2 = 0.058 

ipmax/i < \v\ < Prnax), Mq = 0.8, U = 0.33 {pniax = 1.4, 

dmax = 1.131). Panels (c) and (d) are the results for the initial 
distribution with rji = 0.66 (0 < |p| < Pmax/3), 772 ~ 0.11 

{Prnax/S < \p\ < Pmax) , Mq = 0.9, U = 0.18 {pmax = 1.1, 
dmax = 0.787). 



compare with MD simulations, we plot both the position 
and the velocity distribution functions. As can be seen 
the agreement between the theory and the simulations is 
excellent, without any adjustable parameters. 

It is interesting to compare the predictions of the 
present theory with the approach based on the maximiza- 
tion of the coarse-grained Lynden-Bell (LB) entropy [.18[. 
Within the LB statistics the distribution function is pre- 
dicted to be 
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cos9fsie,p)d0dp = Ms 
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give us a total of 2L + 1 equations to determine M^, and 
the L threshold energies {si} and {xi} values. AU these 
equations have to be solved self-consistently to obtain 
the distribution function in the qSS, Eq. ([7]). In Fig. 
4 we present the solution of these equations for various 
two-level distribution functions satisfying the GVC. To 
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and /3 and {aj} are the Lagrange multipliers used to 
conserve the total energy and the phase space volume 
of each density level. For one-level waterbag distribu- 
tions the LB theory was found to work very well when 
the initial distribution function satisfied the GVC. This, 



however, is no longer the case for the multilevel distri- 
butions. Fig. 4, shows that for such distributions the 
predictions of the LB theory deviate significantly from 
the results of MD simulations. In particular, the LB the- 
ory violates the topological constraint that the stationary 
distribution function must have no microscopic holes in 
the central core region. 



We have studied the dynamics of coUisionless relax- 
ation of systems with LR forces. It was found that if 
the initial phase-space particle distribution has compact 
simply connected support - has no holes - the final sta- 
tionary distribution will also contain a compact region. 
We find that the microscopic holes created by the fila- 
mentation of the initial distribution function are always 
restricted to the "outer" regions of the phase space re- 
sulting in a characteristic core-halo structure. 



For an arbitrary initial multi-level distributions it is 
very difficult to a priori predict the final qSS without 
solving explicitly the full many-body dynamics. An in- 
complete relaxation and a non-ergodic mixing of the dif- 
ferent density levels prevents the use of standard meth- 
ods of statistical mechanics. However, we find that for 
multi-level initial distributions satisfying the GVC it is 
possible to a priori predict the particle distribution in the 
qSS without any adjustable parameters. The challenge 
now is to extend the theory to the initial distributions 
which are not virialized. 
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